#pragma OPENCL EXTENSION cl_khr_fp64 : enable



double xbrz(double t, double t1, double x1, double t2,
            double x2, double r, double xr, double alp, double alk)
{

    if (fabs(r) < 1e-10)
    {
        double ax = x2 - x1;
        double tx = t2 - t1;
        return (x1 + t * ax / tx);
    }
    else
    {
        double alx = alk - alp;
        double tx = t2 - t1;
        double al = alp + t * alx / tx;
        return (xr + fabs(r) * cos(al));
    }
}

double getxbrz(double tbiez, __global double *tn, __global double *xn, __global double *rn, __global double *xrn, __global double *alpn, __global double *alkn)
{
    int i = 1;

    while (tbiez > tn[i] + 1e-10)
    {
        ++i;
    }
    double tloc = tbiez - tn[i - 1];

    return xbrz(tloc, tn[i - 1], xn[i - 1], tn[i], xn[i], rn[i - 1], xrn[i - 1], alpn[i - 1], alkn[i - 1]);

}

double gettx(double t1, double t2, double x, __global double *
             tn, __global double *xn, __global double *rn, __global double *xrn, __global double *alpn,
             __global double *alkn)
{

    double xx, tbiez;
    double tstart = t1;
    double tstop = t2;
    double xstart = getxbrz(tstart, tn, xn, rn, xrn, alpn, alkn);

    while (fabs(tstop - tstart) > 1e-10)
    {
        tbiez = (tstart + tstop) / 2.;
        xx = getxbrz(tbiez, tn, xn, rn, xrn, alpn, alkn);
        if ((x - xx) * (xx - xstart) > 0.)
        {
            tstart = tbiez;
            xstart = xx;
        }
        else
        {
            tstop = tbiez;
        }
    }

    return (tstart + tstop) / 2.;
}

double ybrz0(double t, double t1, double y1, double t2,
             double y2, double r, double yr, double alp, double alk)
{

    if (fabs(r) < 1e-10)
    {
        double ay = y2 - y1;
        double ty = t2 - t1;
        return (y1 + t * ay / ty);
    }
    else
    {
        double aly = alk - alp;
        double ty = t2 - t1;
        double al = alp + t * aly / ty;
        return (yr + fabs(r) * sin(al));
    }
}

double getybrz(double tbiez, __global double *tn, __global double *yn,
               __global double *rn, __global double *yrn, __global double *alpn, __global double *alkn)
{

    int i = 1;
    while (tbiez > tn[i] + 1e-10)
    {
        ++i;
    }

    double tloc = tbiez - tn[i - 1];

    return ybrz0(tloc, tn[i - 1], yn[i - 1], tn[i], yn[i], rn[i - 1], yrn[i - 1], alpn[i - 1], alkn[i - 1]);
}

double gety(double x, double eta, double t1g, double
            t2g, double t1d, double t2d, __global double *tn, __global double *xn,
            __global double *yn, __global double *rn, __global double *xrn, __global double *yrn,
            __global double *alpn, __global double *alkn)
{

    double tbiezg = gettx(t1g, t2g, x, tn, xn, rn, xrn, alpn, alkn);
    double fg = getybrz(tbiezg, tn, yn, rn, yrn, alpn, alkn);
    double tbiezd = gettx(t1d, t2d, x, tn, xn, rn, xrn, alpn, alkn);
    double fd = getybrz(tbiezd, tn, yn, rn, yrn, alpn, alkn);

    return fd + eta * (fg - fd);
}



double getty(double t1, double t2, double y, __global double *
             tn, __global double *yn, __global double *rn, __global double *yrn, __global double *alpn,
             __global double *alkn)
{

    double yy, tbiez;
    double tstart = t1;
    double tstop = t2;
    double ystart = getybrz(tstart, tn, yn, rn, yrn, alpn, alkn);

    while (fabs(tstop - tstart) > 1e-10)
    {
        tbiez = (tstart + tstop) / 2.;
        yy = getybrz(tbiez, tn, yn, rn, yrn, alpn, alkn);
        if ((y - yy) * (yy - ystart) > 0.)
        {
            tstart = tbiez;
            ystart = yy;
        }
        else
        {
            tstop = tbiez;
        }
    }

    return (tstart + tstop) / 2.;
}


double getxbep_(double eta, double xp, double yp,
                double xk, double yk, __global double *xn, __global double *yn,
                __global double *rn, __global double *xrn, __global double *yrn)
{

    double yfu;
    double wsp = xk - xp;
    double etan = xp + eta * wsp;
    double dx, dy, rr, dxb, dyb, cosal, sinal;
    int i = 1;

    if (etan >= xk)
    {
        yfu = yk;
    }
    else if (etan <= xp)
    {
        yfu = yp;
    }
    else
    {

        while (etan < xn[i])
        {
            ++i;
        }
        rr = fabs(rn[i - 1]);
        if (rr < 1e-10)
        {
            dx = xn[i] - xn[i - 1];
            dy = yn[i] - yn[i - 1];
            dxb = etan - xn[i - 1];
            dyb = dy * dxb / dx;
            yfu = yn[i - 1] + dyb;
        }
        else
        {
            sinal = (etan - xrn[i - 1]) / rr;
            cosal = sqrt(1. - sinal * sinal);
            yfu = yrn[i - 1] + rr * cosal;
        }
    }

    return yfu;
}


double getx_(double eta, double ksi, int nn, __global double *
             xn, __global double *yn, __global double *rn, __global double *xrn, __global double *yrn,
             double xcharld, double ycharld, double xcharl,
             double ycharl, int ncharld, int ncharl, double
             xchard, double ychard, double xchar, double ychar,
             int nchard, int nchar)
{


    /*     fl=getxbel(eta,0.d0,xcharld,1.d0,xcharl,yn(2),xn(2), */
    /*    1                                           nn-ncharld+ncharl) */

    //int i__1 = nn - ncharld + ncharl;
    double fl = getxbep_(eta, ycharld, xcharld, ycharl, xcharl, &yn[ncharld - 2],
                         &xn[ncharld - 2], &rn[ncharld - 2], &yrn[ncharld - 2], &xrn[ncharld - 2]);
    //i__1 = nchard - nchar + 1;
    double fp = getxbep_(eta, ychard, xchard, ychar, xchar, &yn[nchar - 2], &xn[nchar - 2],
                         &rn[nchar - 2], &yrn[nchar - 2], &xrn[nchar - 2]);

    return fl + ksi * (fp - fl);
}

double getx_p(double y, double ksi, double t1p,
              double t2p, double t1l, double t2l, __global double *tn,
              __global double *xn, __global double *yn, __global double *rn, __global double *xrn,
              __global double *yrn, __global double *alpn, __global double *alkn)
{
    double tbiezp = getty(t1p, t2p, y, tn, yn, rn, yrn, alpn, alkn);
    double fp = getxbrz(tbiezp, tn, xn, rn, xrn, alpn, alkn);
    double tbiezl = getty(t1l, t2l, y, tn, yn, rn, yrn, alpn, alkn);
    double fl = getxbrz(tbiezl, tn, xn, rn, xrn, alpn, alkn);

    return fl + ksi * (fp - fl);
}


__kernel void podzial(__global double *ypodz, __global double* xpodz, __global double *xn,
                      __global double *yn, __global double *rn, __global double *xrn, __global double *yrn,
                      __global double *alpn, __global double *alkn, __global double *tn, __global double *xnew,
                      int nn, double tcharl, double tchar, double tchard, double tcharld,
                      double xcharld, double ycharld, double xcharl,
                      double ycharl, double xchard,
                      double ychard,double xchar,
                      double ychar, int ncharld, int ncharl, int nchard, int nchar, double tcl, int nw, int mw)
{
    const int i = get_global_id(0);
    const int j = get_global_id(1);

    int j1;
    double xx, yy, xx0, yy0, xsr, ysr;
    int irep;
    
    
    j1 = i*nw + j + 1;
    xx = getx_(ypodz[j], xpodz[i], nn, xn, yn, rn,
               xrn, yrn, xcharld, ycharld, xcharl, ycharl,
               ncharld, ncharl, xchard, ychard, xchar, ychar, nchard, nchar);
    xx0 = 1e10;
    yy = 1e10;
    irep = 0;
    while (fabs(xx - xx0) > 1e-12 && irep < 10)
    {
        ++irep;
        yy0 = yy;
        yy = gety(xx, ypodz[j], tcharl, tchar, tchard, tcharld, tn, xn, yn, rn, xrn, yrn, alpn, alkn);
        xx0 = xx;
        xx = getx_p(yy, xpodz[i], tchar, tchard, tcharld, tcl, tn, xn, yn, rn, xrn, yrn, alpn, alkn);
    }
    if (irep >= 10)
    {
        xsr = xx;
        irep = 0;
        while (fabs(xx - xx0) > 1e-12 && irep < 200)
        {
            ++irep;
            yy0 = yy;
            yy = gety(xsr, ypodz[j], tcharl, tchar, tchard, tcharld, tn, xn, yn, rn, xrn, yrn, alpn, alkn);
            ysr = (yy0 + yy) / 2.;
            xx0 = xx;
            xx = getx_p(ysr, xpodz[i], tchar, tchard, tcharld, tcl, tn, xn, yn, rn, xrn, yrn, alpn, alkn);
            xsr = (xx0 + xx) / 2.;
        }
    }
    xnew[2 * j1 - 2] = xx;
    xnew[2 * j1 - 1] = yy;
}
